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We analyze the structure of the energy landscape of a well-relaxed 1000-atom model of amorphous 
silicon using the activation-relaxation technique (ART nouveau). Generating more than 40,000 
events starting from a single minimum, we find that activated mechanisms are local in nature, that 
they are distributed uniformly throughout the model and that the activation energy is limited by the 
cost of breaking one bond, independently of the complexity of the mechanism. The overall shape of 
the activation-energy-barrier distribution is also insensitive to the exact details of the configuration, 
indicating that well-relaxed configurations see essentially the same environment. These results 
underscore the localized nature of relaxation in this material. 

PACS numbers: 61.43.Dq, 66.30.Hs, 02.70.-C 



I. INTRODUCTION 

The dynamics of many complex materials is dominated 
by activated jumps over energy barriers generally higher 
than ksT. For these systems, the energy-landscape pic- 
ture, which focuses on the topological relation between 
locally stable states, has proven very valuable. Discon- 
nectivity graphs, for example, introduced by Czerminski 
and Elberi and applied extensively by others ; 2 : 3 : 4 : 5 : 6 have 
provided a first classification of the dynamics of complex 
systems based on the structure of their respective energy 
landscape. 

In parallel with these developments, oriented towards 
reconstructing the topology of the energy landscape, 
there has been some efforts in trying to characterize the 
energetics and the nature of the events in clusters t 3 * 4 * 5 ^ 
proteins^*2iS and in amorphous materials. ^SiiLiSii 3 . This 
extensive sampling is essential in order to try to connect 
the properties of the landscape with the dynamics mea- 
sured experimentally, it also serves to build a better un- 
derstanding of the generic properties of various networks: 
low vs. high connectivity, bulk vs. finite-size systems, co- 
valent vs. ionic bonding, etc. 

In this paper, we present an extensive study of the 
structure of the energy landscape of amorphous silicon 
around two minima. Less extensive studies of this mate- 
rial were already presented both by our groupiSiii* 1 ^ and 
Middlcton et a/., 13 using two different approaches. Using 
ART nouveau, we generate more than 42 000 activated 
events around a well-relaxed minimum and analyze the 
properties of the reaction paths connecting this initial 
minimum to a nearby saddle point and new minimum. 
We find that : (1) the activation mechanisms in a-Si are 
local, limiting somewhat the usefulness of the configura- 
tional energy landscape picture, (2) the activation energy 
is essentially limited by the energy required to break a 
single bond, (3) the entropic barrier is almost identical 
for all events, (4) more than 20 % of all events generated 
are bond-switches of the Wooten- Winer- Weaire type, and 
(5) the number of events seems to be of the order 30 to 
60 per atom. 



TABLE I: Parameters of the original Stillinger- Weber poten- 
tial as well as those modified by Vink et al, as described in 
Ref. ITE! We use the latter set in this paper. 
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II. DETAILS OF THE SIMULATION 

We study here the energy landscape around two well- 
relaxed 1000-atom configurations of a-Si. We focus 
mostly on the first one, repeating the simulation on a 
second configuration only to ensure that the results are 
generic and not dependent on some specific feature of the 
configuration. 

The energy model used is a modified version (mSW) 
of the empirical Stillinger- Weber (SW) potentiali 4 * 1 ^ de- 
veloped to ensure that the elastic and structural proper- 
ties of the amorphous phase of silicon are in agreement 
with experiment.- 5 With its original parameters, the SW 
potential describes both the liquid and the crystalline 
phases with good accuracy but fails to reproduce the 
experimental structure of the amorphous phase liS* 1 ^*^ 
Recently, Vink and collaborators proposed a slightly dif- 
ferent set of parameters (see Table HJ, which generates 
the right amorphous structure in addition to providing 
the correct vibrational properties. Because of its em- 
pirical nature and the fact that it is not optimized for 
dynamics, the energy barriers computed with this poten- 
tial must be taken with some care. However, since the 
structural properties of this material are well described 
by mSW, the qualitative features of the energy landscape 
are expected to be correct. 
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A. ART nouveau 

While our previous study of relaxation in a-Si used 
a version of ART that could not identify saddle 
points , 12 i 18 i 19 the results presented here are obtained us- 
ing ART nouveau, the latest version of the activation- 
relaxation technique presented in Ref. 00- Using ideas 
similar to those proposed by Munro and Wales fS^ ART 
nouveau applies a Lanczos scheme to compute directly 
the lowest second derivative of the energy during the acti- 
vation phase and ensures convergence to the saddle point 
to any desired precision. 

Events are generated in the following way: In order 
to leave the harmonic well, one atom is selected at ran- 
dom. This atom and its neighbors, contained in a shell 
of radius 3.5 A, are moved iteratively in a randomly cho- 
sen direction, while keeping the energy, projected in the 
perpendicular directions to a minimum. At each step, a 
Lanczos scheme is used to compute the lowest eigenvalues 
and eigenvectors associated with the curvature of the en- 
ergy landscape. We consider that the configuration has 
left the harmonic well when the lowest eigenvalue falls 
below -50 eV/A 2 . 

The activation process per se starts from this point. 
The configuration is slowly pushed up along the direc- 
tion of lowest curvature until the modulus of the force 
falls below 0.5 eV/A — indicating that the configuration 
has reached a saddle point — or until the lowest eigen- 
value become positive — indicating that the trajectory 
selected is back in the initial harmonic well region and 
must be rejected. After reaching the saddle point, the 
configuration is pushed slightly away from it, and is re- 
laxed into a new local energy minimum, called the '"fi- 
nal minimum" . This event is stored and a new event is 
started from the same initial minimum. 
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FIG. 1: Solid line: radial distribution function (RDF) of the 
1000-atom a-Si model used in the simulation; dashes: exper- 
imental RDF measured by Laaziri el aS~. 



imum in order to ensure that the features of the energy 
landscape were statistically independent of the minimum 
selected. All the quantities analyzed here are the same 
for both minimum, confirming that the properties of the 
landscape are not affected by the details of the topol- 
ogy. In view of the similarities between the properties of 
these two sets of data, most of the discussion will focus 
on Minimum 1. 



B. Properties of the initial configurations 

The two initial configurations used here were prepared 
using ART nouveau and mSW. Starting from a 1000- 
atom randomly packed unit cell with periodic bound- 
ary conditions, ART events were applied until the con- 
figurational energy equilibrated. We used a Metropo- 
lis accept/reject algorithm, based the energy differ- 
ence between consecutive local minima, as described in 
Ref. Il8l2ll The first configuration has an energy per 
atom of -4.000 eV, with 20 five-fold and 26 three-fold 
defects. The radial distribution function is in good agree- 
ment with recent experimental data(Fig. ^) . The model 
is there fore of comparable quality to the models discussed 
in Refs. ll9j22l This configuration is used as the origin of 
all events for the first 42 000-event run. 

The second configuration was prepared by further ap- 
plying ART on the initial configuration, at a Metropolis 
temperature T = 0.25 eV, for a few thousands of events, 
until the average displacement per atom reached 1 A. We 
generated more than 7000 events around this second min- 



C. Search for activated events 



All events presented here are characterized by three 
configurations: the initial minimum (either Minimum 1 
or 2), a first-order saddle point and a second (final) min- 
imum, obtained by relaxing the configurational energy 
from the saddle point, in a direction opposite to that of 
the initial minimum. As discussed in Ref. 0, these events 
are reversible both from the saddle point and the final 
minimum. 

Each event is generated starting from the same initial 
minimum, but with a different random initial direction, 
the rest of the procedure being deterministic. It takes 
about 400 force evaluations to generated a single event. 
About one third of the random launches lead to a trajec- 
tory that brings the configuration back in the harmonic 
basin; these trajectories are simply rejected and a new 
random direction selected. The whole search took about 
4 weeks on a single processor of a Regatta Power4 IBM. 
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III. RESULTS 



Distributions on events 



We first discuss the distributions of events generated 
from Minimum 1. 42 581 events were generated follow- 
ing the procedure discussed in the previous section. The 
distribution of activated (Saddle — initial) and asym- 
metry (-Efi na i — ©initial) energies for all these events are 
given in Fig. [3 The activated energies form a con- 
tinuous spectrum, from to more than 7 eV, with an 
average barrier height of 3.0 eV and a distribution width 
of 1.2 eV. Both the width and average barrier energies 
are much lower than those reported earlier using the pre- 
vious version of ART. 1 ^ This underlines the limit of the 
original ART, which cannot converge directly to a sad- 
dle point. As shown in Ref. II 2t the error on the barrier 
using the initial version of ART could be as high as 1 
eV. As with the method of Munro and Wales^S on the 
contrary, ART nouveau makes it possible to identify the 
transition point to any desired accuracy. Comparing with 
experiment, we find that the average barrier height is in 
overall agreement with measurements of Shin and Atwa- 
te»24 which indicate activation barriers extending from 
0.25 eV to about 2.8. The typical activation energies we 
find correspond also to the isothermal calorimetric data 
of Roorda et at which indicate high- activation barriers. 

Although there are no experimental information to 
compare with, it is useful to examine the distribution 
of the asymmetry energy, i.e., the energy difference be- 
tween the final and initial minima (see the bottom panel 
of Fig. |2J). The average of the distribution is at 1.7 eV, 
with a width of 1.4 eV. There are only a few events with a 
final energy lower than the initial as should be expected, 
since the initial configuration is already very well-relaxed. 
As shown below, the distribution around Minimum 1 is 
essentially identical to that around Minimum 2. This 
results suggests that the overall shape of the configura- 
tional energy landscape is not sensitive to the details of 
the configuration, contrary to what one could think. 

It is also useful to compare our results with those of 
Middleton and Wales obtained on the same system using 
the eigenvector following BFGS approach (EF-BFGS); 13 
ART nouveau is similar in spirit to this method. In 
Fig. |21 we plot the activated energy calculated from the 
initial minimum and the final (©saddle — ©final) for both 
sets of events. The barrier distribution calculated from 
the initial minimum show that the sampling of events dif- 
fers seriously for both methods; while EF-BFGS seems to 
favor strongly events with a barrier below 2 eV, ART se- 
lects events in a more Gaussian way. Since both methods 
follow closely the direction corresponding to the lowest 
eigenvalue to the saddle point, this difference is solely 
due to the algorithm used to leave the harmonic well. 
Remarkably, however, this selection has very little im- 
pact on the barrier distribution calculated from the final 
minimum; in this case, both methods find the same bar- 
rier distribution, which is heavily skewed towards very 
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FIG. 2: Solid lines : Normalized distribution of energies for 
the 42 581 events generated around Minimum 1. Top: dis- 
tribution of the activated energy i? sa ddie — initial- Bottom: 
distribution of the asymmetry energy -E nna i — ^initial- Dotted 
lines : distribution of energies for the 8799 different (unique) 
saddle points (top) and the 6519 different final configurations. 



lower barriers. We discuss the significance of these dis- 
similarities in Section llVl 

While the energy distribution inform us about the bar- 
riers heights, the distribution of Hamming distances be- 
tween configurations (Fig. |2} provide some insight as to 
the rearrangements taking place. As was found previ- 
ously in o-Si and i>-Si02 most saddle points are some- 
where at mid-way between the initial and final stage^Si 
This is somewhat unexpected in the energy landscape 
picture: in a high-dimensional space, any two ran- 
domly selected direction are orthogonal. For truly high- 
dimensional events, therefore, one would expect to find 
little correlation between the displacement at the sad- 
dle point and that at the final minimum. This is not 
the case, however, if events are taking place in a much 
restricted sub-space. The distribution of displacements 
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FIG. 3: Comparison between the distributions of activation 
energy calculated from the final (top) and the initial (bot- 
tom) minimum for the set of events computed by Middleton 
and Wales — - and that generated with ART nouveau from 
minimum 1. 
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FIG. 4: Distribution of the total displacement (Hamming dis- 
tance) for the all events starting from Minimum 1 at the sad- 
dle point (solid line) and at the final minimum (dashed line). 
The distribution of Hamming displacements for unique events, 
as defined in the text, is essentially identical to this one. 



indicates therefore that even though the simulation sys- 
tem is embedded in a 3000-dimensional space, the effec- 
tive sub-space in which each event takes place is much 
smaller and local in nature; one must therefore be cau- 
tious when interpreting the dynamics of a material based 
solely on the energy landscape picture. 
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FIG. 5: Distribution of the number of activated atoms at 
the saddle point (solid line) and the final minimum (dashed 
line) for all events starting at minimum 1. The distribution 
is almost identical when only unique events, as defined in the 
text, are taken into account. 



The local nature of the dynamics is also found plotting 
the number of active atoms participating in each event. 
(Fig.EJ). An active atom is defined as one that has moved 
by more than 0.1 A from its initial state. As discussed 
m Refs. H9H2IL the value of the threshold is chosen to be 
close to the typical atomic displacement due to thermal 
vibration at room temperature. We see that most events 
involve a maximum of 25 atoms at the saddle point, and 
40 at the final minimum. If the threshold is increased to 
0.4 or 0.5 A, the number of active events drops typically 
to between 2 and 4 per event. The relatively narrow dis- 
tribution of active atoms is not an artifact due to some 
bias in the sampling of the cell. Fig. shows the prob- 
ability of each atom to participate into an event. The 
remarkably homogeneous figure confirms that the rear- 
rangements can take place anywhere in the network, with 
a participation probability varying by at most a factor 3. 



B. Properties of the landscape 

As the sampling of events is random, the total data 
set contains redundant events that might bias the analy- 
sis of the landscape. It is therefore useful to look at the 
properties of "unique" events, i.e., to analyze the set of 
single copies of all events. We identify each event based 
on the identity of the participating atoms as well as on 
the energy at the barrier and at the new minimum. In 
order to decrease the impact of the finite precision in 
the convergence at the saddle and the minimum, we de- 
fine the participating atoms as those moving by at least 
<5f threshold = 0.4 and 0.2 A, respectively, at the saddle 
point and the final minimum. Two events are considered 
identical if the active atoms are the same and the energy 
barrier differs by less than 0.2 eV. We verified that the 
precise values of the various thresholds do not affect the 
qualitative results although the precise number of dif- 
ferent events obviously depends on the threshold. The 
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FIG. 6: Probability of participating in a event for each atom 
of the model, computed using on all generated events. Top: 
probability of being an active atom at the saddle point; bot- 
tom: probability of being an active atom at the final mini- 
mum. 
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FIG. 7: Solid line: Number of unique events as a function 
of the number of events sampled. The dashed and the dot- 
ted lines show the number of unique saddle points and final 
minima as a function of the number of events sampled. 



^10° 



3- 
o 

rr 



-2 





2 4 

A E (eV) 



FIG. 8: Top: Log-normal ratio of all saddles points generated 
from minimum 1 over the unique ones only. A histogram for 
both distributions is first constructed, as a function of the 
energy, and the ratio is taken over this histogram. Bottom: 
same but the configurations at the final minimum. The insets 
in both figures shows the same distributions plotted linearly. 



number of unique events, as a function of the number of 
events already sampled, is shown in Fig. We find 6519 
different minima, 8799 different saddles and 11014 differ- 
ent events generated within the 42 000+ sequence. On 
average, each event is therefore visited almost 4 times. 

Using the distribution of unique events, it is possible to 
assess the biases in the ART sampling. Fig. |S] shows the 
ratio of all saddle points or minima generated over the 
list of all unique ones. As was the case for the Lennard- 
Jones clusters, a system with a totally different energy 
landscape, ART sampling seems to select events with an 
exponential bias, exp(— AE/A). 

These results are obtained using a subset of all events 
surrounding Minimum 1. Even after 42000 events, the 
number of different events around minimum is not yet 
converged. It is possible to obtain some very crude es- 
timate of this number. We first set the full distribution 



of saddle points identical to that shown in Fig. [7] Sup- 
posing a random selection of events, with an exponential 
bias on the barrier height such as that of the previous 
paragraph, taken from the distribution of unique events, 
we find that growth of the total of unique events as a 
function of trial can be reproduced with an exponential 
bias exp(— AE/A) with A between 0.40 and 0.60, and a 
total number of events somewhere between 30,000 and 
60,000 or 30 to 60 different saddle points per atom. 



C. Entropy 

Having computed the distribution of activation energy 
barriers around a minimum, it would be possible to con- 
sider associating a time scale with the events. This can 
be done within the framework provided by the transition 
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state theory^. For activated mechanisms, the diffusion 
rate is given by 
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where g is the number of equivalent diffusion paths, a the 
spatial dimension, I the length of the jump, v§ a phonon 
frequency, and AF = AE — TAS, the variation in the 
Helmholtz free energy between the transition state and 
the initial state. 

In order to use the above equation, we need to measure 
the change in entropy from the minimum to the transition 
state. We can do that using the harmonic approximation: 



AS* = k b In 



3N 

n 

i=i 



3JV-1 
i=l 



(2) 



where vf 1 and vf^ are the real phonon frequencies at 
the minimum and at the saddle point respectively. Since 
there are 3N — 1 real frequencies at the transition state, 
we replace the imaginary frequency by a typical phonon 
frequency, Vq. 

Assessing the time scale associated with leaving a given 
minimum remains a very expensive task even within the 
harmonic approximation as it would be necessary to com- 
pute the entropy associated with each of the 30,000 to 
60,000 saddle points. A simpler approach poses that the 
entropic prefactor is independent of the specific event. 

We can ascertain the validity of such assumption by 
computing the distribution of the entropic difference from 
the initial minimum to the saddle point. Due to the 
high cost of diagonalizing a 3000 x 3000 matrix, we have 
applied this procedure to 50 randomly-selected events. 
As seen in Fig. [§| the contribution of the entropy to 
the activation is of the order of 0.0024 eV/K, with a 
variation of about 1.5 x 10~ 4 cV/K. The fluctuations in 
the entropic contributions to the barrier can therefore 
affect the attempt frequency by about a factor 5. To a 
first approximation, it is therefore reasonable to consider 
that the entropy is only a multiplicative constant in the 
dynamics of the system. 



D. Topological classification 

In the previous sub-sections, we have analyzed the 
properties of the energy landscape itself. We now turn 
to the classification of the events generated around a sin- 
gle minimum. The topological classification used here is 
described in detail in Ref. 19. Briefly, it only consid- 
ers changes in the coordination of the atoms, obtained 
through bond-breaking and binding. Analyzing the more 
than 40 000 events generated here, we find more than 
2000 different types of topological events. Of those, we 
discuss below only the events that occurred more than 1 
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FIG. 9: Distribution of the entropy barrier at the saddle point 
evaluated in the harmonic limit, for 50 different events se- 
lected at random. 
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% of the time; about 400 times in our sequence. Only 5 
types of event satisfy this criterion. 

As in previous simulations, the most common event 
we find, representing about 20 % of all generated events, 
is the Wooten- Winer- Weaire bond exchange mechanism 
(or abacbd, in our notation), introduced almost 20 years 
ago in the sillium model. Since then, this mechanism has 
been seen both in crystalline and amorphous materials. 
^*2SiSL22i22^£ Figure ^| shows the energy and Hamming 
distance distribution for this mechanism. Because of the 
frequency of this mechanism, it might not be surprising 
to find that both the energy and Hamming distance dis- 
tributions follow closely those obtained for the whole set 
of events. These distributions underline once more, nev- 
ertheless, the role of strain in determining the activation 
barrier for a given mechanism: the same topological jump 
can lead to a new configuration with an energy varying 
by as much as 6 eV. 

The topological nature of the other 4 dominant events 
is also shown in Fig. ^2 The class abc is associated with 
a bond jumping from on pair of atoms to the next, ab 
represents a bond breaking at the saddle point and re- 
forming. It does not involve any topological change and 
is responsible for the low-energy peak in the asymmetry 
energy distribution. As such it is not very interesting. 
Finally, the two other classes are modifications on the 
three most common types. 



IV. DISCUSSION AND CONCLUSIONS 

The purpose of this study was to characterize in de- 
tail the energy landscape around a single minimum in a 
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FIG. 10: Top: Energy distribution for the 8000 WWW-type 
events. Solid line: distribution at the saddle point; dotted 
line: distribution of the asymmetry energy. Bottom: Distri- 
bution of the Hamming distance to the initial minimum at 
the saddle point (solid line) and the final minimum (dashed 
line) . 
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FIG. 11: The 5 most common types of events, according to 
our topological classification. The full lines represent a bond 
present in the initial configuration and the dotted lines a bond 
present in the final configuration. 



little variation in the entropic barrier in spite of the wide 
spread in the energy part. This is likely related to the 
fact that the overall shape of the energy landscape is 
not finely dependent on the specifics of a configuration. 
Finally, the bias found in ART nouveau indicates that 
it might be possible to contemplate setting up a kinetic 
Monte Carlo simulation on this system. 

It is interesting to note the strong dissimilarities be- 
tween the activation barrier distribution as computed 
from the initial and the final minimum. In particular, 
the barrier measured from the latter are much lower 
that those measured from the former configuration. This 
difference can be explained by the fact that the initial 
configuration selected is very-well relaxed; most low en- 
ergy barrier, associated with a very unstable direction, 
would have been crossed during the relaxation, in pre- 
vious events. These states are likely not contributing 
to diffusion or the relaxation of the network since they 
would rapidly relax back into the initial minimum. The 
physically relevant barrier distribution for defining the 
evolution of the system is therefore that measured from 
the initial minimum. 

In spite of the 40,000+ events generated, the total 
number of events around a minimum remains unknown. 
We can nevertheless estimate that the number of differ- 
ent paths is somewhere between 30 and 60 paths/atom; 
as shown here, this number should be independent of the 
size of the network as all events are local in nature. 

Finally the stability of the distribution from one min- 
imum to another, as well as the narrow distribution of 
entropic barriers suggests that it would be possible to de- 
velop an accelerated algorithm for a-Si, similar to that 
proposed by Hcnkclman and Jonsson for the diffusion of 
Cu on Cu^i, and that of Hernandez-Rojas and Wales for 
LJ glasses. 32 This is made even easier by the exponential 
bias of ART, also seen in LJ clusters. The origin of this 
bias is not understood but it simplifies significantly the 
statistical analysis of an ART sampling. 

At this point, however, the main limitation of this type 
of simulation is the absence of detailed experimental re- 
sults. Although the overall energy scale of the barriers 
is in agreement with experimental numbers, it is difficult 
to assess whether or not the theoretical description is 
correct: there is no information available experimentally 
on the shape of the barrier distribution. It is currently 
impossible to generate a better comparison with experi- 
ment. Further experimental work is clearly needed to ex- 
pand on the current available data and help ensure that 
the advances in numerical results arc on solid grounds. 



well-relaxed model of a-Si. Previous work has examined 
much shorter sequences of events as the configuration re- 
laxed and the atoms diffused. Here, we performed two 
extensive simulations, always starting from the same two 
minima. Our results indicate that events are essentially 
local with a barrier height limited by the cost of break- 
ing one single bond. We also found that there is very 
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